   
REAL(8) FUNCTION func_COMMIT3(x0)

!Performs exercise for forward guidance 
!(transfers are fixed as gN is modified)    
! Value A: Repayment (Baseline) with adjusted current gN
! Value B: Default (Baseline)

!INSTRUCTIONS:
!1. With (na,nb)=(21,61), store baseline solution:
!   - in graphs\base folder: v, q_paid, b_next 
!   - in graphs\base\pol folder: ggn, gh, gp, gtau, gca, gibp
!2. Need to call func beforehand to upload all parameter values
!3. Set parameters below for the specific ayusterity measure
!4. At the end of subroutine IRF3 uncomment only the lines needed to store the relevant results for the exercise conducted
        
USE GLOBAL
USE OMP_LIB

implicit none

REAL(8),INTENT(IN)::x0(6)

INTEGER  i_b, i_a, i_default_global,i0(1),i,j
INTEGER(4)::hours2,mins2,i_b2
REAL(8)::a_initial,b_initial
REAL(8)::stats(6)
REAL(8)::FDATA(nb),DLEFT,DRIGHT,BREAK_GRID(nb),CSCOEF(4,nb),DCSVAL
INTEGER::ILEFT, IRIGHT,NINTV
REAL(8)::new_value
INTEGER::index_bfeas(na),index_opt,ii_b,nbaux,ival
REAL(8)::q_fun,probdef_fun
EXTERNAL::q_fun,probdef_fun

REAL(8)::bcut_grid(nbounds),sizecut_grid(nsizecut),bcutlow,bcuthgh
    
REAL(8)::fvalor,func_gN,xval(6),func
EXTERNAL  func_gN,func

beta =x0(1)
chi  =x0(2)
icost=x0(3)
psi1 =x0(4)
psi2 =x0(5)
wbar =x0(6)

!!Need to upload solution to baseline model
call compute_bgrid

do i=1,nb2
   bgrid2(i) = bmin + (bmax - bmin) * REAL((i-1d0) / (nb2 - 1d0))
end do

do i=1,na2
   agrid2(i) = amin + (amax - amin) * REAL((i-1d0) / (na2 - 1d0))
end do

i0 = LOCATE2(bgrid2,0d0)
i_bzero2 = i0(1)
bgrid2(i_bzero2)=0d0

fvalor=func(x0)
   
index_aglobal=251

DO i=1,1
    DO j=1,1
    
    index_global =LOCATE2(bgrid2,-0.27124d0)

    !current cuts:          size cut=0.03       bcutmin=-10000  bcutmax=10000 
    !blunt cuts:            size cut=0.03       bcutmin=-10000  bcutmax=10000 
    !state-contingent cuts: size cut = 0.04     bcutmin=3.83    bcutmax=5.06
    
    bcutmin = -10000d0 
    bcutmax = 10000d0 
    sizecut = 0.03d0
    
    PRINT*,'bcutmin,bcutmax,sizecut',bcutmin,bcutmax,sizecut

    open (10, FILE='graphs\base\v.txt')             !guess from graphs
    open (16, FILE='graphs\base\q_paid.txt')
    open (17, FILE='graphs\base\b_next.txt')
    open (19, FILE='graphs\base\pol\ggn.txt')       !pol fcns 21 61
    open (20, FILE='graphs\base\pol\gh.txt')
    open (21, FILE='graphs\base\pol\gp.txt')
    open (22, FILE='graphs\base\pol\gtau.txt')
    open (23, FILE='graphs\base\pol\gca.txt')
    open (24, FILE='graphs\base\pol\gibp.txt')

    do i_a = 1,na
    do i_b = 1,nb
        READ(19, *) gn_base(i_a, i_b)   
        READ(20, *) h_base(i_a, i_b)   
        READ(21, *) pricen_base(i_a, i_b)   
        READ(22, *) tau_base(i_a, i_b)   
        READ(23, *) ca_base(i_a, i_b)   
        READ(24, '(I7)') ival
        bp_base(i_a, i_b)=bgrid(ival)   
    ENDDO
    ENDDO

    CLOSE(19)
    CLOSE(20)
    CLOSE(21)
    CLOSE(22)
    CLOSE(23)
    CLOSE(24)

    transf_base = tau_base + ca_base-pricen_base*gn_base
    wage_base = alpha*pricen_base*(h_base**(alpha-1d0))

    do i_b = 1,nb
    do i_a = 1,na
        READ(10, '(F18.10, X, F18.10, X, F15.10)') v_matrix(i_b, i_a), &
                         v0_matrix(i_b, i_a), vaut_matrix(i_b, i_a)
        IF ((na2.EQ.na).AND.(nb2.EQ.nb)) THEN
            gn1(i_a, i_b)=gn_base(i_a, i_b)   
            vcon1(i_a, i_b)=v0_matrix(i_b, i_a)
        ENDIF
        
        READ(16, '(F15.11, X, F15.11)') q_matrix(i_b, i_a), q_matrix_nodef(i_b, i_a)
        READ(17, '(F15.11)') b_next_matrix(i_b, i_a)
        if (vaut_matrix(i_b, i_a) > v0_matrix(i_b, i_a)) then
            default_decision(i_b, i_a) =2
        else
            default_decision(i_b, i_a) =1        
        end if
    end do

    end do


    CLOSE(10)
    CLOSE(16)
    CLOSE(17)

    q_matrixx= q_matrix

    PRINT*, 'Saved Initial Guess Uploaded'

    index_bfeasmax2 = nb+1          !for interim value fcn in period 2
    index_bfeas4 = 0d0
    
    !to compute spline for feasible debt gridpoints
    do i_a = 1,na
    new_value=-1d5
    i_b=0
    index_opt = 0
        
    do while (new_value <= -1d3)
        index_opt = i_b
        IF (i_b.EQ.nb) EXIT
        i_b=i_b+1
        new_value = v0_matrix(i_b,i_a)
    end do
    index_bfeas4(i_a)=index_opt     
    ENDDO

    index_bfeas2=MAXVAL(index_bfeas4)
         
    do i_a = 1,na         
        IF (index_bfeas4(i_a)+1.LE.nb-1) THEN      !need at least to points for splines
            ii_b=index_bfeas4(i_a)+1
            nbaux =nb-ii_b+1
            FDATA(1:nbaux) = v0_matrix(ii_b:nb,i_a)
	        ILEFT  = 0
	        IRIGHT = 0
	        DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(ii_b+1) - bgrid(ii_b))
	        DRIGHT = (FDATA(nbaux)-FDATA(nbaux-1)) / (bgrid(nbaux) - bgrid(nbaux-1))
	        call  DCSDEC (nbaux, bgrid(ii_b:nb), FDATA(1:nbaux), ILEFT, DLEFT, IRIGHT, DRIGHT, break_matrix(i_a, ii_b:nb), coeff_matrix(i_a, :,ii_b:nb))
        ENDIF         
    ENDDO
    
    do i_a = 1,na
        FDATA = q_matrix_nodef(:,i_a)
	    ILEFT  = 0
	    IRIGHT = 0
        DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(2) - bgrid(1))
	    DRIGHT = (FDATA(nb)-FDATA(nb-1)) / (bgrid(nb) - bgrid(nb-1))
	    call  DCSDEC (nb, bgrid, FDATA, ILEFT, DLEFT, IRIGHT, DRIGHT, BREAK_GRID, CSCOEF)
        break_matrix_q(i_a, :) = BREAK_GRID
        coeff_matrix_q(i_a, :,:) = CSCOEF
    end do

    !call val_COMMIT22       !current cut
    !call val_COMMIT23       !next-period cut

    !call val_COMMIT023
    !call val_COMMIT022
    call IRFS3
ENDDO
ENDDO

func_COMMIT3 = ((-stats(1)-22.8d0)/22.8d0)**2d0 +((stats(2)-18.1d0)/18.1d0)**2d0 +((stats(3)-1.1d0)/1.1d0)**2d0 +((stats(4)-1.4d0)/1.4d0)**2d0+((stats(5)-1.11d0)/1.11d0)**2d0+((stats(6)-0.9815d0)/0.9815d0)**2d0

END FUNCTION func_COMMIT3
!***************************************************************************************************************
REAL(8) FUNCTION func_COMMIT3_YT(x0)
!Compute welfare gains, etc. for blunt cuts as function of current yT
        
USE GLOBAL
USE OMP_LIB

implicit none

REAL(8),INTENT(IN)::x0(6)

INTEGER  i_b, i_a, i_default_global,i0(1),i,j,jj
INTEGER(4)::hours2,mins2,i_b2
REAL(8)::a_initial,b_initial
REAL(8)::stats(6)
REAL(8)::FDATA(nb),DLEFT,DRIGHT,BREAK_GRID(nb),CSCOEF(4,nb),DCSVAL
INTEGER::ILEFT, IRIGHT,NINTV
REAL(8)::new_value
INTEGER::index_bfeas(na),index_opt,ii_b,nbaux,ival,index_aux(1)
REAL(8)::q_fun,probdef_fun
EXTERNAL::q_fun,probdef_fun

REAL(8)::bcut_grid(nbounds),sizecut_grid(nsizecut),bcutlow,bcuthgh
    
REAL(8)::fvalor,func_gN,xval(6),func
EXTERNAL  func_gN,func

beta =x0(1)
chi  =x0(2)
icost=x0(3)
psi1 =x0(4)
psi2 =x0(5)
wbar =x0(6)

!!Need to upload solution to baseline model
call compute_bgrid

do i=1,nb2
   bgrid2(i) = bmin + (bmax - bmin) * REAL((i-1d0) / (nb2 - 1d0))
end do

do i=1,na2
   agrid2(i) = amin + (amax - amin) * REAL((i-1d0) / (na2 - 1d0))
end do

i0 = LOCATE2(bgrid2,0d0)
i_bzero2 = i0(1)
bgrid2(i_bzero2)=0d0
      
fvalor=func(x0)

index_aglobal=251


DO i=1,1
    DO j=1,1
    
    index_global =LOCATE2(bgrid2,-0.27124d0)

    !current cuts:          size cut=0.03       bcutmin=-10000  bcutmax=10000 
    !blunt cuts:            size cut=0.03       bcutmin=-10000  bcutmax=10000 
    !state-contingent cuts: size cut = 0.04     bcutmin=3.83    bcutmax=5.06
    
    bcutmin = -10000d0 
    bcutmax = 10000d0 
    sizecut = 0.03d0   
    
    PRINT*,'bcutmin,bcutmax,sizecut',bcutmin,bcutmax,sizecut

    open (10, FILE='graphs\base\v.txt')             !guess from graphs
    open (16, FILE='graphs\base\q_paid.txt')
    open (17, FILE='graphs\base\b_next.txt')
    open (19, FILE='graphs\base\pol\ggn.txt')       
    open (20, FILE='graphs\base\pol\gh.txt')
    open (21, FILE='graphs\base\pol\gp.txt')
    open (22, FILE='graphs\base\pol\gtau.txt')
    open (23, FILE='graphs\base\pol\gca.txt')
    open (24, FILE='graphs\base\pol\gibp.txt')

    do i_a = 1,na
    do i_b = 1,nb
        READ(19, *) gn_base(i_a, i_b)   
        READ(20, *) h_base(i_a, i_b)   
        READ(21, *) pricen_base(i_a, i_b)   
        READ(22, *) tau_base(i_a, i_b)   
        READ(23, *) ca_base(i_a, i_b)   
        READ(24, '(I7)') ival
        bp_base(i_a, i_b)=bgrid(ival)   
    ENDDO
    ENDDO

    CLOSE(19)
    CLOSE(20)
    CLOSE(21)
    CLOSE(22)
    CLOSE(23)
    CLOSE(24)

    transf_base = tau_base + ca_base-pricen_base*gn_base
    wage_base = alpha*pricen_base*(h_base**(alpha-1d0))

    do i_b = 1,nb
    do i_a = 1,na
        READ(10, '(F18.10, X, F18.10, X, F15.10)') v_matrix(i_b, i_a), &
                         v0_matrix(i_b, i_a), vaut_matrix(i_b, i_a)
        IF ((na2.EQ.na).AND.(nb2.EQ.nb)) THEN
            gn1(i_a, i_b)=gn_base(i_a, i_b)   
            vcon1(i_a, i_b)=v0_matrix(i_b, i_a)
        ENDIF
        
        READ(16, '(F15.11, X, F15.11)') q_matrix(i_b, i_a), q_matrix_nodef(i_b, i_a)
        READ(17, '(F15.11)') b_next_matrix(i_b, i_a)
        if (vaut_matrix(i_b, i_a) > v0_matrix(i_b, i_a)) then
            default_decision(i_b, i_a) =2
        else
            default_decision(i_b, i_a) =1        
        end if
    end do

    end do


    CLOSE(10)
    CLOSE(16)
    CLOSE(17)

    q_matrixx= q_matrix

    PRINT*, 'Saved Initial Guess Uploaded'

    index_bfeasmax2 = nb+1          !for interim value fcn in period 2
    index_bfeas4 = 0d0
    
    !to compute spline for feasible debt gridpoints
    do i_a = 1,na
    new_value=-1d5
    i_b=0
    index_opt = 0
        
    do while (new_value <= -1d3)
        index_opt = i_b
        IF (i_b.EQ.nb) EXIT
        i_b=i_b+1
        new_value = v0_matrix(i_b,i_a)
    end do
    index_bfeas4(i_a)=index_opt     
    ENDDO

    index_bfeas2=MAXVAL(index_bfeas4)
         
    do i_a = 1,na         
        IF (index_bfeas4(i_a)+1.LE.nb-1) THEN      !need at least to points for splines
            ii_b=index_bfeas4(i_a)+1
            nbaux =nb-ii_b+1
            FDATA(1:nbaux) = v0_matrix(ii_b:nb,i_a)
	        ILEFT  = 0
	        IRIGHT = 0
	        DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(ii_b+1) - bgrid(ii_b))
	        DRIGHT = (FDATA(nbaux)-FDATA(nbaux-1)) / (bgrid(nbaux) - bgrid(nbaux-1))
	        call  DCSDEC (nbaux, bgrid(ii_b:nb), FDATA(1:nbaux), ILEFT, DLEFT, IRIGHT, DRIGHT, break_matrix(i_a, ii_b:nb), coeff_matrix(i_a, :,ii_b:nb))
        ENDIF         
    ENDDO
    
    do i_a = 1,na
        FDATA = q_matrix_nodef(:,i_a)
	    ILEFT  = 0
	    IRIGHT = 0
        DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(2) - bgrid(1))
	    DRIGHT = (FDATA(nb)-FDATA(nb-1)) / (bgrid(nb) - bgrid(nb-1))
	    call  DCSDEC (nb, bgrid, FDATA, ILEFT, DLEFT, IRIGHT, DRIGHT, BREAK_GRID, CSCOEF)
        break_matrix_q(i_a, :) = BREAK_GRID
        coeff_matrix_q(i_a, :,:) = CSCOEF
    end do
    
        !call val_COMMIT22       !current cut
        call val_COMMIT23       !next-period cut

        !as function of current yT
        DO jj=1,na
            index_aux=MINLOC((agrid2-agrid(jj))**2d0)
            index_aglobal=index_aux(1)
            PRINT*,'a jj',jj,agrid2(index_aglobal),agrid(jj)

            call val_COMMIT023
            
            call IRFS3
        ENDDO    
ENDDO
ENDDO

func_COMMIT3_YT = ((-stats(1)-22.8d0)/22.8d0)**2d0 +((stats(2)-18.1d0)/18.1d0)**2d0 +((stats(3)-1.1d0)/1.1d0)**2d0 +((stats(4)-1.4d0)/1.4d0)**2d0+((stats(5)-1.11d0)/1.11d0)**2d0+((stats(6)-0.9815d0)/0.9815d0)**2d0

END FUNCTION func_COMMIT3_YT
!***************************************************************************************************************
REAL(8) FUNCTION func_COMMIT3_debt(x0)
!Compute welfare gains, etc. for blunt cuts as function of current debt
        
USE GLOBAL
USE OMP_LIB

implicit none

REAL(8),INTENT(IN)::x0(6)

INTEGER  i_b, i_a, i_default_global,i0(1),i,j,jj
INTEGER(4)::hours2,mins2,i_b2
REAL(8)::a_initial,b_initial
REAL(8)::stats(6)
REAL(8)::FDATA(nb),DLEFT,DRIGHT,BREAK_GRID(nb),CSCOEF(4,nb),DCSVAL
INTEGER::ILEFT, IRIGHT,NINTV
REAL(8)::new_value
INTEGER::index_bfeas(na),index_opt,ii_b,nbaux,ival,index_aux(1)
REAL(8)::q_fun,probdef_fun
EXTERNAL::q_fun,probdef_fun

REAL(8)::bcut_grid(nbounds),sizecut_grid(nsizecut),bcutlow,bcuthgh
    
REAL(8)::fvalor,func_gN,func,xval(6)
EXTERNAL  func_gN,func

beta =x0(1)
chi  =x0(2)
icost=x0(3)
psi1 =x0(4)
psi2 =x0(5)
wbar =x0(6)

!!Need to upload solution to baseline model
call compute_bgrid

do i=1,nb2
   bgrid2(i) = bmin + (bmax - bmin) * REAL((i-1d0) / (nb2 - 1d0))
end do

do i=1,na2
   agrid2(i) = amin + (amax - amin) * REAL((i-1d0) / (na2 - 1d0))
end do

i0 = LOCATE2(bgrid2,0d0)
i_bzero2 = i0(1)
bgrid2(i_bzero2)=0d0
      
fvalor=func(x0)

index_aglobal=251

DO i=1,1 
    DO j=1,1
    
   !fvalor=func_gN(x0)
    index_global = 3172 +j 
    
    !current cuts:          size cut=0.03       bcutmin=-10000  bcutmax=10000 
    !blunt cuts:            size cut=0.03       bcutmin=-10000  bcutmax=10000 
    !state-contingent cuts: size cut = 0.04     bcutmin=3.83    bcutmax=5.06
    
    bcutmin = -10000d0 
    bcutmax = 100000d0 
    sizecut = 0.03d0   
    
    PRINT*,'bcutmin,bcutmax,sizecut',bcutmin,bcutmax,sizecut

    open (10, FILE='graphs\base\v.txt')             !guess from graphs
    open (16, FILE='graphs\base\q_paid.txt')
    open (17, FILE='graphs\base\b_next.txt')
    open (19, FILE='graphs\base\pol\ggn.txt')      
    open (20, FILE='graphs\base\pol\gh.txt')
    open (21, FILE='graphs\base\pol\gp.txt')
    open (22, FILE='graphs\base\pol\gtau.txt')
    open (23, FILE='graphs\base\pol\gca.txt')
    open (24, FILE='graphs\base\pol\gibp.txt')


    do i_a = 1,na
    do i_b = 1,nb
        READ(19, *) gn_base(i_a, i_b)   
        READ(20, *) h_base(i_a, i_b)   
        READ(21, *) pricen_base(i_a, i_b)   
        READ(22, *) tau_base(i_a, i_b)   
        READ(23, *) ca_base(i_a, i_b)   
        READ(24, '(I7)') ival
        bp_base(i_a, i_b)=bgrid(ival)   
    ENDDO
    ENDDO

    CLOSE(19)
    CLOSE(20)
    CLOSE(21)
    CLOSE(22)
    CLOSE(23)
    CLOSE(24)

    transf_base = tau_base + ca_base-pricen_base*gn_base
    wage_base = alpha*pricen_base*(h_base**(alpha-1d0))

    do i_b = 1,nb
    do i_a = 1,na
        READ(10, '(F18.10, X, F18.10, X, F15.10)') v_matrix(i_b, i_a), &
                         v0_matrix(i_b, i_a), vaut_matrix(i_b, i_a)
        IF ((na2.EQ.na).AND.(nb2.EQ.nb)) THEN
            gn1(i_a, i_b)=gn_base(i_a, i_b)   
            vcon1(i_a, i_b)=v0_matrix(i_b, i_a)
        ENDIF
        
        READ(16, '(F15.11, X, F15.11)') q_matrix(i_b, i_a), q_matrix_nodef(i_b, i_a)
        READ(17, '(F15.11)') b_next_matrix(i_b, i_a)
        if (vaut_matrix(i_b, i_a) > v0_matrix(i_b, i_a)) then
            default_decision(i_b, i_a) =2
        else
            default_decision(i_b, i_a) =1        
        end if
    end do
    end do


    CLOSE(10)
    CLOSE(16)
    CLOSE(17)

    q_matrixx= q_matrix

    PRINT*, 'Saved Initial Guess Uploaded'

    index_bfeasmax2 = nb+1          !for interim value fcn in period 2
    index_bfeas4 = 0d0
    
    !to compute spline for feasible debt gridpoints
    do i_a = 1,na
    new_value=-1d5
    i_b=0
    index_opt = 0
        
    do while (new_value <= -1d3)
        index_opt = i_b
        IF (i_b.EQ.nb) EXIT
        i_b=i_b+1
        new_value = v0_matrix(i_b,i_a)
    end do
    index_bfeas4(i_a)=index_opt     
    ENDDO

    index_bfeas2=MAXVAL(index_bfeas4)
         
    do i_a = 1,na         
        IF (index_bfeas4(i_a)+1.LE.nb-1) THEN      !need at least to points for splines
            ii_b=index_bfeas4(i_a)+1
            nbaux =nb-ii_b+1
            FDATA(1:nbaux) = v0_matrix(ii_b:nb,i_a)
	        ILEFT  = 0
	        IRIGHT = 0
	        DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(ii_b+1) - bgrid(ii_b))
	        DRIGHT = (FDATA(nbaux)-FDATA(nbaux-1)) / (bgrid(nbaux) - bgrid(nbaux-1))
	        call  DCSDEC (nbaux, bgrid(ii_b:nb), FDATA(1:nbaux), ILEFT, DLEFT, IRIGHT, DRIGHT, break_matrix(i_a, ii_b:nb), coeff_matrix(i_a, :,ii_b:nb))
        ENDIF         
    ENDDO
    
    do i_a = 1,na
        FDATA = q_matrix_nodef(:,i_a)
	    ILEFT  = 0
	    IRIGHT = 0
        DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(2) - bgrid(1))
	    DRIGHT = (FDATA(nb)-FDATA(nb-1)) / (bgrid(nb) - bgrid(nb-1))
	    call  DCSDEC (nb, bgrid, FDATA, ILEFT, DLEFT, IRIGHT, DRIGHT, BREAK_GRID, CSCOEF)
        break_matrix_q(i_a, :) = BREAK_GRID
        coeff_matrix_q(i_a, :,:) = CSCOEF
    end do

    call val_COMMIT23      !next-period cut

    DO jj=8,nb
        index_aux=MINLOC((bgrid2-bgrid(jj))**2d0)
        index_global=index_aux(1)
        PRINT*,'b jj',jj,bgrid2(index_global),bgrid(jj)

        call val_COMMIT023
        call IRFS3
    ENDDO    
ENDDO
ENDDO

func_COMMIT3_debt = ((-stats(1)-22.8d0)/22.8d0)**2d0 +((stats(2)-18.1d0)/18.1d0)**2d0 +((stats(3)-1.1d0)/1.1d0)**2d0 +((stats(4)-1.4d0)/1.4d0)**2d0+((stats(5)-1.11d0)/1.11d0)**2d0+((stats(6)-0.9815d0)/0.9815d0)**2d0

END FUNCTION func_COMMIT3_debt
!***************************************************************************************************************
REAL(8) FUNCTION func_COMMIT3_SIZECUT(x0)
!Compute welfare gains for blunt cuts and state-contingent cuts depending on size of cut        

USE GLOBAL
USE OMP_LIB

implicit none

REAL(8),INTENT(IN)::x0(6)

INTEGER  i_b, i_a, i_default_global,i0(1),i,j
INTEGER(4)::hours2,mins2,i_b2
REAL(8)::a_initial,b_initial
REAL(8)::stats(6)
REAL(8)::FDATA(nb),DLEFT,DRIGHT,BREAK_GRID(nb),CSCOEF(4,nb),DCSVAL
INTEGER::ILEFT, IRIGHT,NINTV
REAL(8)::new_value
INTEGER::index_bfeas(na),index_opt,ii_b,nbaux,ival
REAL(8)::q_fun,probdef_fun
EXTERNAL::q_fun,probdef_fun

REAL(8)::bcut_grid(nbounds),sizecut_grid(nsizecut),bcutlow,bcuthgh
    
REAL(8)::fvalor,func_gN,func,xval(6)
EXTERNAL  func_gN,func

beta =x0(1)
chi  =x0(2)
icost=x0(3)
psi1 =x0(4)
psi2 =x0(5)
wbar =x0(6)

!!Need to upload solution to baseline model
call compute_bgrid

do i=1,nb2
   bgrid2(i) = bmin + (bmax - bmin) * REAL((i-1d0) / (nb2 - 1d0))
end do

do i=1,na2
   agrid2(i) = amin + (amax - amin) * REAL((i-1d0) / (na2 - 1d0))
end do

i0 = LOCATE2(bgrid2,0d0)
i_bzero2 = i0(1)
bgrid2(i_bzero2)=0d0

fvalor=func(x0)

!size of spending cuts:
sizecut_grid(1)=0.01d0
sizecut_grid(2)=0.02d0
sizecut_grid(3)=0.03d0
sizecut_grid(4)=0.04d0
sizecut_grid(5)=0.05d0
sizecut_grid(6)=0.06d0

index_aglobal=251

DO i=1,1 
    DO j=1,6
    
   !fvalor=func_gN(x0)
    index_global =LOCATE2(bgrid2,-0.27124d0) 
    
    !Set bcutmin and bcutmax depending on blunt cuts vs. state-contingent cuts
    !blunt cuts:            bcutmin=-10000  bcutmax=10000
    !state-contingent cuts: bcutmin=3.83    bcutmax=5.06
    
    bcutmin = 3.83d0 
    bcutmax = 5.06d0 
    sizecut = sizecut_grid(j)
    
    PRINT*,'bcutmin,bcutmax,sizecut',bcutmin,bcutmax,sizecut

    open (10, FILE='graphs\base\v.txt')             !guess from graphs
    open (16, FILE='graphs\base\q_paid.txt')
    open (17, FILE='graphs\base\b_next.txt')
    open (19, FILE='graphs\base\pol\ggn.txt')       
    open (20, FILE='graphs\base\pol\gh.txt')
    open (21, FILE='graphs\base\pol\gp.txt')
    open (22, FILE='graphs\base\pol\gtau.txt')
    open (23, FILE='graphs\base\pol\gca.txt')
    open (24, FILE='graphs\base\pol\gibp.txt')

    do i_a = 1,na
    do i_b = 1,nb
        READ(19, *) gn_base(i_a, i_b)   
        READ(20, *) h_base(i_a, i_b)   
        READ(21, *) pricen_base(i_a, i_b)   
        READ(22, *) tau_base(i_a, i_b)   
        READ(23, *) ca_base(i_a, i_b)   
        READ(24, '(I7)') ival
        bp_base(i_a, i_b)=bgrid(ival)   
    ENDDO
    ENDDO

    CLOSE(19)
    CLOSE(20)
    CLOSE(21)
    CLOSE(22)
    CLOSE(23)
    CLOSE(24)

    transf_base = tau_base + ca_base-pricen_base*gn_base
    wage_base = alpha*pricen_base*(h_base**(alpha-1d0))

    do i_b = 1,nb
    do i_a = 1,na
        READ(10, '(F18.10, X, F18.10, X, F15.10)') v_matrix(i_b, i_a), &
                         v0_matrix(i_b, i_a), vaut_matrix(i_b, i_a)
        IF ((na2.EQ.na).AND.(nb2.EQ.nb)) THEN
            gn1(i_a, i_b)=gn_base(i_a, i_b)   
            vcon1(i_a, i_b)=v0_matrix(i_b, i_a)
        ENDIF
        
        READ(16, '(F15.11, X, F15.11)') q_matrix(i_b, i_a), q_matrix_nodef(i_b, i_a)
        READ(17, '(F15.11)') b_next_matrix(i_b, i_a)
        if (vaut_matrix(i_b, i_a) > v0_matrix(i_b, i_a)) then
            default_decision(i_b, i_a) =2
        else
            default_decision(i_b, i_a) =1        
        end if
    end do

    end do


    CLOSE(10)
    CLOSE(16)
    CLOSE(17)

    q_matrixx= q_matrix

    PRINT*, 'Saved Initial Guess Uploaded'

    index_bfeasmax2 = nb+1          !for interim value fcn in period 2
    index_bfeas4 = 0d0
    
    !to compute spline for feasible debt gridpoints
    do i_a = 1,na
    new_value=-1d5
    i_b=0
    index_opt = 0
        
    do while (new_value <= -1d3)
        index_opt = i_b
        IF (i_b.EQ.nb) EXIT
        i_b=i_b+1
        new_value = v0_matrix(i_b,i_a)
    end do
    index_bfeas4(i_a)=index_opt     
    ENDDO

    index_bfeas2=MAXVAL(index_bfeas4)
         
    do i_a = 1,na         
        IF (index_bfeas4(i_a)+1.LE.nb-1) THEN      !need at least to points for splines
            ii_b=index_bfeas4(i_a)+1
            nbaux =nb-ii_b+1
            FDATA(1:nbaux) = v0_matrix(ii_b:nb,i_a)
	        ILEFT  = 0
	        IRIGHT = 0
	        DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(ii_b+1) - bgrid(ii_b))
	        DRIGHT = (FDATA(nbaux)-FDATA(nbaux-1)) / (bgrid(nbaux) - bgrid(nbaux-1))
	        call  DCSDEC (nbaux, bgrid(ii_b:nb), FDATA(1:nbaux), ILEFT, DLEFT, IRIGHT, DRIGHT, break_matrix(i_a, ii_b:nb), coeff_matrix(i_a, :,ii_b:nb))
        ENDIF         
    ENDDO
    
    do i_a = 1,na
        FDATA = q_matrix_nodef(:,i_a)
	    ILEFT  = 0
	    IRIGHT = 0
        DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(2) - bgrid(1))
	    DRIGHT = (FDATA(nb)-FDATA(nb-1)) / (bgrid(nb) - bgrid(nb-1))
	    call  DCSDEC (nb, bgrid, FDATA, ILEFT, DLEFT, IRIGHT, DRIGHT, BREAK_GRID, CSCOEF)
        break_matrix_q(i_a, :) = BREAK_GRID
        coeff_matrix_q(i_a, :,:) = CSCOEF
    end do
    
    call val_COMMIT23      !next-period cut
    call val_COMMIT023
    call IRFS3
ENDDO
ENDDO

func_COMMIT3_SIZECUT = ((-stats(1)-22.8d0)/22.8d0)**2d0 +((stats(2)-18.1d0)/18.1d0)**2d0 +((stats(3)-1.1d0)/1.1d0)**2d0 +((stats(4)-1.4d0)/1.4d0)**2d0+((stats(5)-1.11d0)/1.11d0)**2d0+((stats(6)-0.9815d0)/0.9815d0)**2d0

END FUNCTION func_COMMIT3_SIZECUT
!**********************************************************************************
REAL(8) FUNCTION func_COMMIT3_LOWERBOUND(x0)
!Compute welfare gains for state-contingent cuts depending on lower bound for income
USE GLOBAL
USE OMP_LIB

implicit none

REAL(8),INTENT(IN)::x0(6)

INTEGER  i_b, i_a, i_default_global,i0(1),i,j
INTEGER(4)::hours2,mins2,i_b2
REAL(8)::a_initial,b_initial
REAL(8)::stats(6)
REAL(8)::FDATA(nb),DLEFT,DRIGHT,BREAK_GRID(nb),CSCOEF(4,nb),DCSVAL
INTEGER::ILEFT, IRIGHT,NINTV
REAL(8)::new_value
INTEGER::index_bfeas(na),index_opt,ii_b,nbaux,ival
REAL(8)::q_fun,probdef_fun
EXTERNAL::q_fun,probdef_fun

REAL(8)::bcut_grid(nbounds),sizecut_grid(nsizecut),bcutlow,bcuthgh
    
REAL(8)::fvalor,func_gN,func,xval(6)
EXTERNAL  func_gN,func

beta =x0(1)
chi  =x0(2)
icost=x0(3)
psi1 =x0(4)
psi2 =x0(5)
wbar =x0(6)

!!Need to upload solution to baseline model
call compute_bgrid

do i=1,nb2
   bgrid2(i) = bmin + (bmax - bmin) * REAL((i-1d0) / (nb2 - 1d0))
end do

do i=1,na2
   agrid2(i) = amin + (amax - amin) * REAL((i-1d0) / (na2 - 1d0))
end do

i0 = LOCATE2(bgrid2,0d0)
i_bzero2 = i0(1)
bgrid2(i_bzero2)=0d0

fvalor=func(x0)
      
bcut_grid(1) = 3.5d0
bcut_grid(2)= 3.55d0
bcut_grid(3) = 3.6d0
bcut_grid(4)= 3.65d0
bcut_grid(5) = 3.7d0
bcut_grid(6)= 3.75d0
bcut_grid(7) = 3.8d0
bcut_grid(8)= 3.85d0
bcut_grid(9) = 3.9d0
bcut_grid(10) = 3.95d0
bcut_grid(11) = 4d0
bcut_grid(12) = 4.05d0
bcut_grid(13) = 4.1d0

index_aglobal=251

DO i=1,1 
    DO j=1,13

    index_global =LOCATE2(bgrid2,-0.27124d0) 
    
    !Set bcutmin and bcutmax depending on blunt cuts vs. state-contingent cuts
    !blunt cuts:            bcutmin=-10000  bcutmax=10000
    !state-contingent cuts: bcutmin=3.83    bcutmax=5.06

    bcutmin = bcut_grid(j)  
    bcutmax = 5.06d0  
    sizecut = 0.04d0
    
    PRINT*,'bcutmin,bcutmax,sizecut',bcutmin,bcutmax,sizecut

    open (10, FILE='graphs\base\v.txt')             !guess from graphs
    open (16, FILE='graphs\base\q_paid.txt')
    open (17, FILE='graphs\base\b_next.txt')
    open (19, FILE='graphs\base\pol\ggn.txt')     
    open (20, FILE='graphs\base\pol\gh.txt')
    open (21, FILE='graphs\base\pol\gp.txt')
    open (22, FILE='graphs\base\pol\gtau.txt')
    open (23, FILE='graphs\base\pol\gca.txt')
    open (24, FILE='graphs\base\pol\gibp.txt')

    do i_a = 1,na
    do i_b = 1,nb
        READ(19, *) gn_base(i_a, i_b)   
        READ(20, *) h_base(i_a, i_b)   
        READ(21, *) pricen_base(i_a, i_b)   
        READ(22, *) tau_base(i_a, i_b)   
        READ(23, *) ca_base(i_a, i_b)   
        READ(24, '(I7)') ival
        bp_base(i_a, i_b)=bgrid(ival)   
    ENDDO
    ENDDO

    CLOSE(19)
    CLOSE(20)
    CLOSE(21)
    CLOSE(22)
    CLOSE(23)
    CLOSE(24)

    transf_base = tau_base + ca_base-pricen_base*gn_base
    wage_base = alpha*pricen_base*(h_base**(alpha-1d0))

    do i_b = 1,nb
    do i_a = 1,na
        READ(10, '(F18.10, X, F18.10, X, F15.10)') v_matrix(i_b, i_a), &
                         v0_matrix(i_b, i_a), vaut_matrix(i_b, i_a)
        IF ((na2.EQ.na).AND.(nb2.EQ.nb)) THEN
            gn1(i_a, i_b)=gn_base(i_a, i_b)   
            vcon1(i_a, i_b)=v0_matrix(i_b, i_a)
        ENDIF
        
        READ(16, '(F15.11, X, F15.11)') q_matrix(i_b, i_a), q_matrix_nodef(i_b, i_a)
        READ(17, '(F15.11)') b_next_matrix(i_b, i_a)
        if (vaut_matrix(i_b, i_a) > v0_matrix(i_b, i_a)) then
            default_decision(i_b, i_a) =2
        else
            default_decision(i_b, i_a) =1        
        end if
    end do

    end do


    CLOSE(10)
    CLOSE(16)
    CLOSE(17)

    q_matrixx= q_matrix

    PRINT*, 'Saved Initial Guess Uploaded'

    index_bfeasmax2 = nb+1          !for interim value fcn in period 2
    index_bfeas4 = 0d0
    
    !to compute spline for feasible debt gridpoints
    do i_a = 1,na
    new_value=-1d5
    i_b=0
    index_opt = 0
        
    do while (new_value <= -1d3)
        index_opt = i_b
        IF (i_b.EQ.nb) EXIT
        i_b=i_b+1
        new_value = v0_matrix(i_b,i_a)
    end do
    index_bfeas4(i_a)=index_opt     
    ENDDO

    index_bfeas2=MAXVAL(index_bfeas4)
         
    do i_a = 1,na         
        IF (index_bfeas4(i_a)+1.LE.nb-1) THEN      !need at least to points for splines
            ii_b=index_bfeas4(i_a)+1
            nbaux =nb-ii_b+1
            FDATA(1:nbaux) = v0_matrix(ii_b:nb,i_a)
	        ILEFT  = 0
	        IRIGHT = 0
	        DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(ii_b+1) - bgrid(ii_b))
	        DRIGHT = (FDATA(nbaux)-FDATA(nbaux-1)) / (bgrid(nbaux) - bgrid(nbaux-1))
	        call  DCSDEC (nbaux, bgrid(ii_b:nb), FDATA(1:nbaux), ILEFT, DLEFT, IRIGHT, DRIGHT, break_matrix(i_a, ii_b:nb), coeff_matrix(i_a, :,ii_b:nb))
        ENDIF         
    ENDDO
    
    do i_a = 1,na
        FDATA = q_matrix_nodef(:,i_a)
	    ILEFT  = 0
	    IRIGHT = 0
        DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(2) - bgrid(1))
	    DRIGHT = (FDATA(nb)-FDATA(nb-1)) / (bgrid(nb) - bgrid(nb-1))
	    call  DCSDEC (nb, bgrid, FDATA, ILEFT, DLEFT, IRIGHT, DRIGHT, BREAK_GRID, CSCOEF)
        break_matrix_q(i_a, :) = BREAK_GRID
        coeff_matrix_q(i_a, :,:) = CSCOEF
    end do

    
    call val_COMMIT23      !next-period cut
    call val_COMMIT023
    call IRFS3
ENDDO
ENDDO

func_COMMIT3_LOWERBOUND = ((-stats(1)-22.8d0)/22.8d0)**2d0 +((stats(2)-18.1d0)/18.1d0)**2d0 +((stats(3)-1.1d0)/1.1d0)**2d0 +((stats(4)-1.4d0)/1.4d0)**2d0+((stats(5)-1.11d0)/1.11d0)**2d0+((stats(6)-0.9815d0)/0.9815d0)**2d0

END FUNCTION func_COMMIT3_LOWERBOUND
!**********************************************************************************
REAL(8) FUNCTION func_COMMIT3_UPPERBOUND(x0)
!Compute welfare gains for state-contingent cuts depending on lower bound for income

USE GLOBAL
USE OMP_LIB

implicit none

REAL(8),INTENT(IN)::x0(6)

INTEGER  i_b, i_a, i_default_global,i0(1),i,j
INTEGER(4)::hours2,mins2,i_b2
REAL(8)::a_initial,b_initial
REAL(8)::stats(6)
REAL(8)::FDATA(nb),DLEFT,DRIGHT,BREAK_GRID(nb),CSCOEF(4,nb),DCSVAL
INTEGER::ILEFT, IRIGHT,NINTV
REAL(8)::new_value
INTEGER::index_bfeas(na),index_opt,ii_b,nbaux,ival
REAL(8)::q_fun,probdef_fun
EXTERNAL::q_fun,probdef_fun

REAL(8)::bcut_grid(17),sizecut_grid(nsizecut),bcutlow,bcuthgh,func
    
REAL(8)::fvalor,func_gN,xval(6)
EXTERNAL  func_gN,func

beta =x0(1)
chi  =x0(2)
icost=x0(3)
psi1 =x0(4)
psi2 =x0(5)
wbar =x0(6)

!!Need to upload solution to baseline model
call compute_bgrid

do i=1,nb2
   bgrid2(i) = bmin + (bmax - bmin) * REAL((i-1d0) / (nb2 - 1d0))
end do

do i=1,na2
   agrid2(i) = amin + (amax - amin) * REAL((i-1d0) / (na2 - 1d0))
end do

i0 = LOCATE2(bgrid2,0d0)
i_bzero2 = i0(1)
bgrid2(i_bzero2)=0d0
      

fvalor=func(x0)

bcut_grid(1)  = 4.5d0
bcut_grid(2)  = 4.55d0
bcut_grid(3)  = 4.6d0
bcut_grid(4)  = 4.65d0
bcut_grid(5)  = 4.7d0
bcut_grid(6)  = 4.75d0
bcut_grid(7)  = 4.8d0
bcut_grid(8)  = 4.85d0
bcut_grid(9)  = 4.9d0
bcut_grid(10) = 4.95d0
bcut_grid(11) = 5.0d0
bcut_grid(12) = 5.05d0
bcut_grid(13) = 5.1d0
bcut_grid(14) = 5.15d0
bcut_grid(15) = 5.2d0
bcut_grid(16) = 5.25d0
bcut_grid(17) = 5.3d0

index_aglobal=251

DO i=1,1 
    DO j=1,17
    
    index_global =LOCATE2(bgrid2,-0.27124d0) 
    
    !Set bcutmin and bcutmax depending on blunt cuts vs. state-contingent cuts
    !blunt cuts:            bcutmin=-10000  bcutmax=10000
    !state-contingent cuts: bcutmin=3.83    bcutmax=5.06

    bcutmin = 3.83d0    
    bcutmax = bcut_grid(j)  
    sizecut = 0.04d0

    open (10, FILE='graphs\base\v.txt')             !guess from graphs
    open (16, FILE='graphs\base\q_paid.txt')
    open (17, FILE='graphs\base\b_next.txt')
    open (19, FILE='graphs\base\pol\ggn.txt')      
    open (20, FILE='graphs\base\pol\gh.txt')
    open (21, FILE='graphs\base\pol\gp.txt')
    open (22, FILE='graphs\base\pol\gtau.txt')
    open (23, FILE='graphs\base\pol\gca.txt')
    open (24, FILE='graphs\base\pol\gibp.txt')

    do i_a = 1,na
    do i_b = 1,nb
        READ(19, *) gn_base(i_a, i_b)   
        READ(20, *) h_base(i_a, i_b)   
        READ(21, *) pricen_base(i_a, i_b)   
        READ(22, *) tau_base(i_a, i_b)   
        READ(23, *) ca_base(i_a, i_b)   
        READ(24, '(I7)') ival
        bp_base(i_a, i_b)=bgrid(ival)   
    ENDDO
    ENDDO

    CLOSE(19)
    CLOSE(20)
    CLOSE(21)
    CLOSE(22)
    CLOSE(23)
    CLOSE(24)

    transf_base = tau_base + ca_base-pricen_base*gn_base
    wage_base = alpha*pricen_base*(h_base**(alpha-1d0))

    do i_b = 1,nb
    do i_a = 1,na
        READ(10, '(F18.10, X, F18.10, X, F15.10)') v_matrix(i_b, i_a), &
                         v0_matrix(i_b, i_a), vaut_matrix(i_b, i_a)
        IF ((na2.EQ.na).AND.(nb2.EQ.nb)) THEN
            gn1(i_a, i_b)=gn_base(i_a, i_b)   
            vcon1(i_a, i_b)=v0_matrix(i_b, i_a)
        ENDIF
        
        READ(16, '(F15.11, X, F15.11)') q_matrix(i_b, i_a), q_matrix_nodef(i_b, i_a)
        READ(17, '(F15.11)') b_next_matrix(i_b, i_a)
        if (vaut_matrix(i_b, i_a) > v0_matrix(i_b, i_a)) then
            default_decision(i_b, i_a) =2
        else
            default_decision(i_b, i_a) =1        
        end if
    end do

    end do


    CLOSE(10)
    CLOSE(16)
    CLOSE(17)

    q_matrixx= q_matrix

    PRINT*, 'Saved Initial Guess Uploaded'

    index_bfeasmax2 = nb+1          !for interim value fcn in period 2
    index_bfeas4 = 0d0
    
    !to compute spline for feasible debt gridpoints
    do i_a = 1,na
    new_value=-1d5
    i_b=0
    index_opt = 0
        
    do while (new_value <= -1d3)
        index_opt = i_b
        IF (i_b.EQ.nb) EXIT
        i_b=i_b+1
        new_value = v0_matrix(i_b,i_a)
    end do
    index_bfeas4(i_a)=index_opt     
    ENDDO

    index_bfeas2=MAXVAL(index_bfeas4)
         
    do i_a = 1,na         
        IF (index_bfeas4(i_a)+1.LE.nb-1) THEN      !need at least to points for splines
            ii_b=index_bfeas4(i_a)+1
            nbaux =nb-ii_b+1
            FDATA(1:nbaux) = v0_matrix(ii_b:nb,i_a)
	        ILEFT  = 0
	        IRIGHT = 0
	        DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(ii_b+1) - bgrid(ii_b))
	        DRIGHT = (FDATA(nbaux)-FDATA(nbaux-1)) / (bgrid(nbaux) - bgrid(nbaux-1))
	        call  DCSDEC (nbaux, bgrid(ii_b:nb), FDATA(1:nbaux), ILEFT, DLEFT, IRIGHT, DRIGHT, break_matrix(i_a, ii_b:nb), coeff_matrix(i_a, :,ii_b:nb))
        ENDIF         
    ENDDO
    
    do i_a = 1,na
        FDATA = q_matrix_nodef(:,i_a)
	    ILEFT  = 0
	    IRIGHT = 0
        DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(2) - bgrid(1))
	    DRIGHT = (FDATA(nb)-FDATA(nb-1)) / (bgrid(nb) - bgrid(nb-1))
	    call  DCSDEC (nb, bgrid, FDATA, ILEFT, DLEFT, IRIGHT, DRIGHT, BREAK_GRID, CSCOEF)
        break_matrix_q(i_a, :) = BREAK_GRID
        coeff_matrix_q(i_a, :,:) = CSCOEF
    end do
    
    call val_COMMIT23      !next-period cut
    call val_COMMIT023
    call IRFS3
ENDDO
ENDDO

func_COMMIT3_UPPERBOUND = ((-stats(1)-22.8d0)/22.8d0)**2d0 +((stats(2)-18.1d0)/18.1d0)**2d0 +((stats(3)-1.1d0)/1.1d0)**2d0 +((stats(4)-1.4d0)/1.4d0)**2d0+((stats(5)-1.11d0)/1.11d0)**2d0+((stats(6)-0.9815d0)/0.9815d0)**2d0

END FUNCTION func_COMMIT3_UPPERBOUND
!**********************************************************************************
subroutine val_COMMIT23

USE GLOBAL
USE OMP_LIB

IMPLICIT NONE

INTEGER :: d

REAL(8)::a_initial,b_initial
INTEGER::i_default_global
REAL(8)::a_valor,b_valor,b_valor2,v_valor,q_fun,b_next,b0_next
REAL(8)::b_next1(na2),ytotal_base(na,nb),spread_base(na,nb)

REAL(8)::q,interpolate,b_next0,interpolate3
REAL(8)::FDATA(nb),DLEFT,DRIGHT,BREAK_GRID(nb),CSCOEF(4,nb),DCSVAL

INTEGER i_b, i_a, i_g,i, j, ILEFT, IRIGHT,NINTV,DNORIN
REAL(8)::x0,b_fun,bmin_feas
REAL(8)::new_value,vaut_matrix_2(nb,na),v0_matrix_2(nb,na)
INTEGER::index_bfeas(na),index_opt,ii_b,nbaux,try2,try
CHARACTER(len=29)::filename
REAL(8),DIMENSION(nb,na)::q_matrix_nodef_2,v_matrix_2,q_scheme_2,probdef_matrix_2
REAL(8)::objective_fun,objective_fun_gN,probdef_fun
EXTERNAL::objective_fun,objective_fun_gN,probdef_fun
REAL(8)::gn_optimal,dif1,dif2,gnmin1,gnmax1
REAL(8)::xq,xb,xh,blevel
REAL(8)::haut_matrix_2(na),gnaut_matrix_2(na),ctaut_matrix_2(na)
REAL(8)::h_matrix_2(nb,na),gn_matrix_2(nb,na),ct_matrix_2(nb,na)
 
DO i_a=1,na
    ytotal_base(i_a,:) = DEXP(agrid(i_a))+pricen_base(i_a,:)*h_base(i_a,:)**alpha 
    spread_base(i_a,:) = (1d0+q_matrix_nodef(:,i_a)**(-1d0)-lambda)/(1d0+r)-1d0
ENDDO
            
!to compute spline for feasible debt gridpoints
do i_a = 1,na
    new_value=-1d5
    i_b=0
    index_opt = 0
        
    do while (new_value <= -1d3)
        index_opt = i_b
        IF (i_b.EQ.nb) EXIT
        i_b=i_b+1
        new_value = v0_matrix(i_b,i_a)
    end do
    index_bfeas4(i_a)=index_opt     
ENDDO

index_bfeas2=MAXVAL(index_bfeas4)
    
do i_a = 1,na         
    IF (index_bfeas4(i_a)+1.LE.nb-1) THEN      !need at least to points for splines
        ii_b=index_bfeas4(i_a)+1
        nbaux =nb-ii_b+1
        FDATA(1:nbaux) = v0_matrix(ii_b:nb,i_a)
	    ILEFT  = 0
	    IRIGHT = 0
	    DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(ii_b+1) - bgrid(ii_b))
	    DRIGHT = (FDATA(nbaux)-FDATA(nbaux-1)) / (bgrid(nbaux) - bgrid(nbaux-1))
	    call  DCSDEC (nbaux, bgrid(ii_b:nb), FDATA(1:nbaux), ILEFT, DLEFT, IRIGHT, DRIGHT, break_matrix(i_a, ii_b:nb), coeff_matrix(i_a, :,ii_b:nb))
    ENDIF         
ENDDO
         
do i_a = 1,na
    FDATA = q_matrix_nodef(:,i_a)
	ILEFT  = 0
	IRIGHT = 0
    DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(2) - bgrid(1))
	DRIGHT = (FDATA(nb)-FDATA(nb-1)) / (bgrid(nb) - bgrid(nb-1))
	call  DCSDEC (nb, bgrid, FDATA, ILEFT, DLEFT, IRIGHT, DRIGHT, BREAK_GRID, CSCOEF)
    break_matrix_q(i_a, :) = BREAK_GRID
    coeff_matrix_q(i_a, :,:) = CSCOEF
end do
        
PRINT*, 'SOLVING FOR COMMITMENT in Period 2'

!To solve model
    IF (psi1.LE.500d0) THEN
    !! Compute the value of repayment
    LOOP_A1: do i_a = 1,na
                  a_initial   = agrid(i_a)
                  b_initial   = bgrid(1)
                  
                  i_default_global=2 !Country defaults
                  b_valor = 0d0

                  !choose gN optimally in autarky
                  v_valor= -objective_fun(b_valor,a_initial,b_initial,i_default_global)
                  
                  b_next1(i_a)   = b_valor
                  vaut_matrix_2(:, i_a) = v_valor
                  haut_matrix_2(i_a)  = x_global
                  gnaut_matrix_2(i_a) = gn_global
                  ctaut_matrix_2(i_a) = ctrad_global
         end do LOOP_A1
    ELSE
            b_next1 = 0d0
            vaut_matrix_2 = -1d10
    ENDIF
        
    h_matrix_2  = 0d0
    
    LOOP_A2: do i_a = 1,na
        blevel=1d2
    LOOP_B:  do i_b = 1,nb 
        a_initial = agrid(i_a)
        b_initial = bgrid(i_b)
        i_default_global=1 !Country does not default
        try=0
        
661         bp_global = b_next_matrix(i_b, i_a)
        
       !in repayment, set gN equal to constant_gN and transfers to equilibrium values
       !if ytotal belongs to (bcutmin,bcutmax)
       IF ((ytotal_base(i_a, i_b).GT.bcutmin).AND.((ytotal_base(i_a,i_b).LE.bcutmax))) THEN         
            constant_gn = (1d0-sizecut)*gn_base(i_a, i_b) 
            transf_global = transf_base(i_a, i_b) 
            call optimize_gN2(b_next0, v_valor,a_initial,b_initial,i_default_global)
       !in repayment, set gN and transfers to equilibrium values
       !if ytotal does not belong to (bcutmin,bcutmax)
       ELSE
           constant_gn = gn_base(i_a, i_b) 
           transf_global = transf_base(i_a, i_b) 
           call optimize(b_next0, v_valor,a_initial,b_initial,i_default_global)
       ENDIF
        
        IF (constant_gn.LT.1d-6) THEN
            v_valor=-1d6
            v0_matrix_2(i_b, i_a) = -1d6
            x0     = 0d0
            b_next0 = 0d0
            h_matrix_2(i_b,i_a)  = 0d0
            gn_matrix_2(i_b,i_a) = 0d0
            ct_matrix_2(i_b,i_a) = 0d0
            
            q_matrix_nodef_2(i_b, i_a) = 0d0            
        ELSE             
            v0_matrix_2(i_b, i_a) = v_valor                  
            x0 = x_global
            h_matrix_2(i_b,i_a)  = x_global
            gn_matrix_2(i_b,i_a) = gn_global
            ct_matrix_2(i_b,i_a) = ctrad_global
        ENDIF

        q_matrix_nodef_2(i_b, i_a) = q_fun(b_next0, a_initial)
        !probdef_matrix_2(i_b, i_a) = probdef_fun(b_next0, a_initial)
        b_next_matrix_2(i_b, i_a) = b_next0   
            
        if (vaut_matrix_2(i_b, i_a) > v0_matrix_2(i_b, i_a)) THEN 
            default_decision_2(i_b,i_a) = 2
            v_matrix_2(i_b, i_a) = vaut_matrix_2(i_b, i_a)
        else                                          
            default_decision_2(i_b, i_a) = 1
            v_matrix_2(i_b, i_a) = v0_matrix_2(i_b, i_a)
        END if

        q_scheme_2(i_b, i_a) =  q_fun(b_initial,a_initial)
        q_matrix_2(i_b, i_a) = q_fun(b_next_matrix_2(i_b,i_a),a_initial)  

    PRINT*,i_b,i_a,v0_matrix_2(i_b, i_a),v_matrix_2(i_b, i_a) 

    end do LOOP_B
    end do LOOP_A2   
    
    !save solution for period 2 in subfolder graphs/commit/period2
    open (10, FILE='graphs/commit/period2/v_matrix.txt',STATUS='replace')
!    open (15, FILE='graphs/commit/period2/probdef_matrix.txt', STATUS='replace')
    open (16, FILE='graphs/commit/period2/q_paid_matrix.txt', STATUS='replace')
    open (17, FILE='graphs/commit/period2/gn_matrix.txt', STATUS='replace')
    open (18, FILE='graphs/commit/period2/h_matrix.txt',STATUS='replace')
    open (19, FILE='graphs/commit/period2/b_next_matrix.txt',STATUS='replace')
    open (20, FILE='graphs/commit/period2/ct_matrix.txt',STATUS='replace')
    open (21, FILE='graphs/commit/period2/def_matrix.txt',STATUS='replace')
      
    do i_b = 1,nb
        do i_a = 1,na
            WRITE(10, '(F15.8, X, F15.8, X, F15.8)') MAX(-99d0,v_matrix_2(i_b, i_a)), &
            MAX(-99d0,v0_matrix_2(i_b, i_a)), MAX(-99d0,vaut_matrix_2(i_b, i_a))
            !WRITE(15, '(F15.11)') probdef_matrix_2(i_b, i_a)
            WRITE(16, '(F15.11, X, F15.11)') q_matrix_2(i_b, i_a), q_matrix_nodef_2(i_b, i_a)
            WRITE(17, '(F15.11)') gn_matrix_2(i_b, i_a)
            WRITE(18, '(F15.11)') h_matrix_2(i_b, i_a)
            WRITE(19, '(F15.11)') b_next_matrix_2(i_b, i_a)
            WRITE(20, '(F15.11)') ct_matrix_2(i_b, i_a)
            WRITE(21, '(I7)') default_decision_2(i_b, i_a)
        end do
    end do

    CLOSE(10)
    !CLOSE(15)
    CLOSE(16)
    CLOSE(17)
    CLOSE(18)
    CLOSE(19)
    CLOSE(20)
    CLOSE(21)

RETURN


PRINT*, 'VALUES FOR COMMITMENT in Period 2 solved'

end subroutine val_COMMIT23
!**********************************************************************************
!optimal policy fcns in period 2
subroutine val_COMMIT22

USE GLOBAL
USE OMP_LIB

IMPLICIT NONE

INTEGER :: d

REAL(8)::a_initial,b_initial
INTEGER::i_default_global
REAL(8)::a_valor,b_valor,b_valor2,v_valor,q_fun,b_next,b0_next
REAL(8)::b_next1(na2),ytotal_base(na,nb),spread_base(na,nb)

REAL(8)::q,interpolate,b_next0,interpolate3
REAL(8)::FDATA(nb),DLEFT,DRIGHT,BREAK_GRID(nb),CSCOEF(4,nb),DCSVAL

INTEGER i_b, i_a, i_g,i, j, ILEFT, IRIGHT,NINTV,DNORIN
REAL(8)::x0,b_fun,bmin_feas
REAL(8)::new_value,vaut_matrix_2(nb,na),v0_matrix_2(nb,na)
INTEGER::index_bfeas(na),index_opt,ii_b,nbaux,try2,try
CHARACTER(len=29)::filename
REAL(8),DIMENSION(nb,na)::q_matrix_nodef_2,v_matrix_2,q_scheme_2,probdef_matrix_2
REAL(8)::objective_fun,objective_fun_gN,probdef_fun
EXTERNAL::objective_fun,objective_fun_gN,probdef_fun
REAL(8)::gn_optimal,dif1,dif2,gnmin1,gnmax1
REAL(8)::xq,xb,xh
REAL(8)::haut_matrix_2(na),gnaut_matrix_2(na),ctaut_matrix_2(na)
REAL(8)::h_matrix_2(nb,na),gn_matrix_2(nb,na),ct_matrix_2(nb,na)
 
DO i_a=1,na
    ytotal_base(i_a,:) = DEXP(agrid(i_a))+pricen_base(i_a,:)*h_base(i_a,:)**alpha 
    spread_base(i_a,:) = (1d0+q_matrix_nodef(:,i_a)**(-1d0)-lambda)/(1d0+r)-1d0
ENDDO
            
!to compute spline for feasible debt gridpoints
do i_a = 1,na
    new_value=-1d5
    i_b=0
    index_opt = 0
        
    do while (new_value <= -1d3)
        index_opt = i_b
        IF (i_b.EQ.nb) EXIT
        i_b=i_b+1
        new_value = v0_matrix(i_b,i_a)
    end do
    index_bfeas4(i_a)=index_opt     
ENDDO

index_bfeas2=MAXVAL(index_bfeas4)
    
do i_a = 1,na         
    IF (index_bfeas4(i_a)+1.LE.nb-1) THEN      !need at least to points for splines
        ii_b=index_bfeas4(i_a)+1
        nbaux =nb-ii_b+1
        FDATA(1:nbaux) = v0_matrix(ii_b:nb,i_a)
	    ILEFT  = 0
	    IRIGHT = 0
	    DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(ii_b+1) - bgrid(ii_b))
	    DRIGHT = (FDATA(nbaux)-FDATA(nbaux-1)) / (bgrid(nbaux) - bgrid(nbaux-1))
	    call  DCSDEC (nbaux, bgrid(ii_b:nb), FDATA(1:nbaux), ILEFT, DLEFT, IRIGHT, DRIGHT, break_matrix(i_a, ii_b:nb), coeff_matrix(i_a, :,ii_b:nb))
    ENDIF         
ENDDO
         
do i_a = 1,na
    FDATA = q_matrix_nodef(:,i_a)
	ILEFT  = 0
	IRIGHT = 0
    DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(2) - bgrid(1))
	DRIGHT = (FDATA(nb)-FDATA(nb-1)) / (bgrid(nb) - bgrid(nb-1))
	call  DCSDEC (nb, bgrid, FDATA, ILEFT, DLEFT, IRIGHT, DRIGHT, BREAK_GRID, CSCOEF)
    break_matrix_q(i_a, :) = BREAK_GRID
    coeff_matrix_q(i_a, :,:) = CSCOEF
end do
        
PRINT*, 'SOLVING FOR COMMITMENT in Period 2'

!To solve model
    IF (psi1.LE.500d0) THEN
    !! Compute the value of repayment
    LOOP_A1: do i_a = 1,na
                  a_initial   = agrid(i_a)
                  b_initial   = bgrid(1)
                  
                  i_default_global=2 !Country defaults
                  b_valor = 0d0

                  !set gN and transfers optimally in autarky
                  v_valor= -objective_fun(b_valor,a_initial,b_initial,i_default_global)
                  
                  b_next1(i_a)   = b_valor
                  vaut_matrix_2(:, i_a) = v_valor
                  haut_matrix_2(i_a)  = x_global
                  gnaut_matrix_2(i_a) = gn_global
                  ctaut_matrix_2(i_a) = ctrad_global
         end do LOOP_A1
    ELSE
            b_next1 = 0d0
            vaut_matrix_2 = -1d10
    ENDIF
        
    h_matrix_2  = 0d0
    
    LOOP_A2: do i_a = 1,na
    LOOP_B:  do i_b = 1,nb 
        a_initial = agrid(i_a)
        b_initial = bgrid(i_b)
        i_default_global=1 !Country does not default
        try=0
        
661         bp_global = b_next_matrix(i_b, i_a)
        
        !set gN and transfers to equilibrium values in repayment in second period for first-period spending cuts
        constant_gn = gn_base(i_a, i_b) 
        transf_global = transf_base(i_a, i_b) 
        call optimize_gN2(b_next0, v_valor,a_initial,b_initial,i_default_global)

        PRINT*,i_a,i_b,v_valor,b_next0

        IF (constant_gn.LT.1d-6) THEN
            v_valor=-1d6
            v0_matrix_2(i_b, i_a) = -1d6
            x0     = 0d0
            b_next0 = 0d0
            h_matrix_2(i_b,i_a)  = 0d0
            gn_matrix_2(i_b,i_a) = 0d0
            ct_matrix_2(i_b,i_a) = 0d0
            !probdef_matrix_2(i_b, i_a) = 1d0
        ELSE             
            v0_matrix_2(i_b, i_a) = v_valor             
            x0 = x_global
            h_matrix_2(i_b,i_a)  = x_global
            gn_matrix_2(i_b,i_a) = gn_global
            ct_matrix_2(i_b,i_a) = ctrad_global
        ENDIF

        q_matrix_nodef_2(i_b, i_a) = q_fun(b_next0, a_initial)
        !probdef_matrix_2(i_b, i_a) = probdef_fun(b_next0, a_initial)
        b_next_matrix_2(i_b, i_a) = b_next0  
            
        if (vaut_matrix_2(i_b, i_a) > v0_matrix_2(i_b, i_a)) THEN !IT IS OPTIMAL TO DEFAULT
            default_decision_2(i_b,i_a) = 2
            v_matrix_2(i_b, i_a) = vaut_matrix_2(i_b, i_a)
        else                                           
            default_decision_2(i_b, i_a) = 1
            v_matrix_2(i_b, i_a) = v0_matrix_2(i_b, i_a)
        END if

        q_scheme_2(i_b, i_a) =  q_fun(b_initial,a_initial)
        q_matrix_2(i_b, i_a) = q_fun(b_next_matrix_2(i_b,i_a),a_initial)  
        
    end do LOOP_B
    end do LOOP_A2   
    
    !save solution for period 2 in subfolder graphs/commit/period2
    open (10, FILE='graphs/commit/period2/v_matrix.txt',STATUS='replace')
!    open (15, FILE='graphs/commit/period2/probdef_matrix.txt', STATUS='replace')
    open (16, FILE='graphs/commit/period2/q_paid_matrix.txt', STATUS='replace')
    open (17, FILE='graphs/commit/period2/gn_matrix.txt', STATUS='replace')
    open (18, FILE='graphs/commit/period2/h_matrix.txt',STATUS='replace')
    open (19, FILE='graphs/commit/period2/b_next_matrix.txt',STATUS='replace')
    open (20, FILE='graphs/commit/period2/ct_matrix.txt',STATUS='replace')
    open (21, FILE='graphs/commit/period2/def_matrix.txt',STATUS='replace')
      
    do i_b = 1,nb
        do i_a = 1,na
            WRITE(10, '(F15.8, X, F15.8, X, F15.8)') MAX(-99d0,v_matrix_2(i_b, i_a)), &
            MAX(-99d0,v0_matrix_2(i_b, i_a)), MAX(-99d0,vaut_matrix_2(i_b, i_a))
            !WRITE(15, '(F15.11)') probdef_matrix_2(i_b, i_a)
            WRITE(16, '(F15.11, X, F15.11)') q_matrix_2(i_b, i_a), q_matrix_nodef_2(i_b, i_a)
            WRITE(17, '(F15.11)') gn_matrix_2(i_b, i_a)
            WRITE(18, '(F15.11)') h_matrix_2(i_b, i_a)
            WRITE(19, '(F15.11)') b_next_matrix_2(i_b, i_a)
            WRITE(20, '(F15.11)') ct_matrix_2(i_b, i_a)
            WRITE(21, '(I7)') default_decision_2(i_b, i_a)
        end do
    end do

    CLOSE(10)
    !CLOSE(15)
    CLOSE(16)
    CLOSE(17)
    CLOSE(18)
    CLOSE(19)
    CLOSE(20)
    CLOSE(21)

RETURN

PRINT*, 'VALUES FOR COMMITMENT in Period 2 solved'

end subroutine val_COMMIT22
!**********************************************************************************
subroutine val_COMMIT023

USE GLOBAL
!USE OMP_LIB

IMPLICIT NONE

INTEGER :: d

REAL(8)::a_initial,b_initial
INTEGER::i_default_global
REAL(8)::a_valor,b_valor,b_valor2,v_valor,v_valor2,q_fun,b_next,b0_next
REAL(8)::b_next1(na2) 

REAL(8)::q,interpolate,b_next0
REAL(8)::FDATA(nb),DLEFT,DRIGHT,BREAK_GRID(nb),CSCOEF(4,nb),DCSVAL

INTEGER i_b, i_a, i_g,i, j, ILEFT, IRIGHT,NINTV,DNORIN
REAL(8)::x0,x1(na),b_fun,bmin_feas,new_value
INTEGER::index_bfeas(na),index_opt,ii_b,iii_b,nbaux,load_results
INTEGER::i_gmax,i_gmax0(1),i_a_index,i_b_index
CHARACTER(len=29)::filename
REAL(8),DIMENSION(nb,na)::q_matrix_nodef_2,v_matrix_2,q_scheme_2,v0_matrix_2,vaut_matrix_2,probdef_matrix_2
REAL(8),DIMENSION(nb2,na2)::v_matrix_20,q_scheme_20,q_matrix_nodef_20
INTEGER,DIMENSION(nb,na)::index_2
REAL(8)::objective_fun,objective_fun_2,interpolate3,b00,probdef_fun
EXTERNAL::objective_fun,objective_fun_2,probdef_fun

load_results = 1
PRINT*,'PERIOD 1: UPLOADING VALUE FCNS AND PRICES FOR PERIOD 2'

!upload solution for period 2
IF (load_results.EQ.1) THEN
    open (10, FILE='graphs/commit/period2/v_matrix.txt',STATUS='OLD')
    !open (15, FILE='graphs/commit/period2/probdef_matrix.txt', STATUS='OLD')
    open (16, FILE='graphs/commit/period2/q_paid_matrix.txt', STATUS='OLD')
      
    do i_b = 1,nb
    do i_a = 1,na                
        READ(10, '(F15.8, X, F15.8, X, F15.8)') v_matrix_2(i_b, i_a), v0_matrix_2(i_b, i_a), vaut_matrix_2(i_b, i_a)
        !READ(15, '(F15.11)') probdef_matrix_2(i_b, i_a)
        READ(16, '(F15.11, X, F15.11)') q_matrix_2(i_b, i_a), q_matrix_nodef_2(i_b, i_a)
    end do
    end do
    CLOSE(10)
    !CLOSE(15)
    CLOSE(16)
ENDIF
       
!to compute value fcns and prices correctly, set
v0_matrix(:,:)  =v0_matrix_2(:,:)
vaut_matrix(:,:)=vaut_matrix_2(:,:)
q_matrix_nodef(:,:) = q_matrix_nodef_2(:,:)
        
!to compute spline for feasible debt gridpoints
do i_a = 1,na
    new_value=-1d5
    i_b=0
    index_opt = 0
        
    do while (new_value <= -1d3)
        index_opt = i_b
        IF (i_b.EQ.nb) EXIT
        i_b=i_b+1
        new_value = v0_matrix(i_b,i_a)
    end do
    index_bfeas4(i_a)=index_opt     
ENDDO

index_bfeas2=MAXVAL(index_bfeas4)

do i_a = 1,na         
    IF (index_bfeas4(i_a)+1.LE.nb-1) THEN      !need at least to points for splines
        ii_b=index_bfeas4(i_a)+1
        nbaux =nb-ii_b+1
        FDATA(1:nbaux) = v0_matrix(ii_b:nb,i_a)
	    ILEFT  = 0
	    IRIGHT = 0
	    DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(ii_b+1) - bgrid(ii_b))
	    DRIGHT = (FDATA(nbaux)-FDATA(nbaux-1)) / (bgrid(nbaux) - bgrid(nbaux-1))
	    call  DCSDEC (nbaux, bgrid(ii_b:nb), FDATA(1:nbaux), ILEFT, DLEFT, IRIGHT, DRIGHT, break_matrix(i_a, ii_b:nb), coeff_matrix(i_a, :,ii_b:nb))
    ENDIF         
ENDDO
         
do i_a = 1,na
    FDATA = q_matrix_nodef(:,i_a)
	ILEFT  = 0
	IRIGHT = 0
    DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(2) - bgrid(1))
	DRIGHT = (FDATA(nb)-FDATA(nb-1)) / (bgrid(nb) - bgrid(nb-1))
	call  DCSDEC (nb, bgrid, FDATA, ILEFT, DLEFT, IRIGHT, DRIGHT, BREAK_GRID, CSCOEF)
    break_matrix_q(i_a, :) = BREAK_GRID
    coeff_matrix_q(i_a, :,:) = CSCOEF
end do


IF (psi1.LE.500d0) THEN
    !! Compute the value of repayment
    LOOP_A1: do i_a = 1,na2
        a_initial   = agrid2(i_a)
        b_initial   = bgrid2(1)
                  
        i_default_global=2 !Country defaults
        b_valor = 0d0

        !choose gN and transfers optimally in autarky
        v_valor= -objective_fun(b_valor,a_initial,b_initial,i_default_global)
        haut1_20(i_a)  = x_global
        gnaut1_20(i_a) = gn_global
        ctaut1_20(i_a) = ctrad_global
        b_next1(i_a)   = b_valor
        vaut1_20(i_a)  = v_valor
    end do LOOP_A1
ELSE
    haut1_20  = 0d0
    gnaut1_20 = 0d0
    ctaut1_20 = 0d0
    b_next1    = 0d0
    vaut1_20  = -1d10
ENDIF

i_a_index =index_aglobal
i_b_index = index_global 

    !!!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(a_initial,b_initial,i_default_global,v_valor,b_next0)
    !!!$OMP DO COLLAPSE(2)
    LOOP_A2: do i_a = i_a_index,i_a_index 
    LOOP_B:  do i_b = i_b_index,i_b_index
                  a_initial = agrid2(i_a)
                  b_initial = bgrid2(i_b)

                  i_default_global=1  !Country does not default
                  
                  !choose gN and transfers optimally in period 1 in repayment for period-2 spending cuts
                  call optimize(b_next0, v_valor,a_initial,b_initial,i_default_global)
                  PRINT*,i_a,i_b,a_initial,v_valor
                  vcon1_20(i_b, i_a) = v_valor               

                  x0 = x_global

                gn1_20(i_b, i_a) = gn_global   
                h1_20(i_b, i_a)  = x0                     
                ct1_20(i_b, i_a) = ctrad_global   
                b1_20(i_b, i_a) = b_next0  
                        
                  if (vaut1_20(i_a) > vcon1_20(i_b, i_a)) THEN 
                     def1_20(i_b,i_a) = 0
                     v1_20(i_b, i_a) = vaut1_20(i_a)
                  else
                    def1_20(i_b, i_a) = 1
                    v1_20(i_b, i_a) = vcon1_20(i_b, i_a)
                  END if

                  q1_20(i_b, i_a) = q_fun(b1_20(i_b, i_a),a_initial) 

    end do LOOP_B
    end do LOOP_A2
    !!!$OMP END DO
    !!!$OMP END PARALLEL
    
    !save solution for first period in subfolder graphs/commit/period1        
    open (10, FILE='graphs/commit/period1/v.txt',STATUS='UNKNOWN')
    open (11, FILE='graphs/commit/period1/default.txt',STATUS='UNKNOWN')
    open (12, FILE='graphs/commit/period1/bb_next.txt',STATUS='UNKNOWN')
    open (13, FILE='graphs/commit/period1/b_next.txt',STATUS='UNKNOWN')
    open (16, FILE='graphs/commit/period1/q_paid.txt', STATUS='UNKNOWN')
    open (114, FILE='graphs/commit/period1/x.txt',STATUS='UNKNOWN')
    open (115, FILE='graphs/commit/period1/gn.txt',STATUS='UNKNOWN')
    open (116, FILE='graphs/commit/period1/ctt.txt',STATUS='UNKNOWN')
    open (117, FILE='graphs/commit/period1/ign.txt',STATUS='UNKNOWN')
    
    do i_b = 1,nb2
    do i_a = 1,na2
        WRITE(10, '(F15.8, X, F15.8, X, F15.8)') MAX(-99d0,v1_20(i_b, i_a)), &
        MAX(-99d0,vcon1_20(i_b, i_a)), MAX(-99d0,vaut1_20(i_a))
        WRITE(11, '(I7)') def1_20(i_b, i_a)
        bp1_20(i_b,i_a)    = LOCATE2(bgrid2,b1_20(i_b, i_a))
        
        WRITE(12, '(F15.11)') b1_20(i_b, i_a)
        WRITE(13, '(I7)') bp1_20(i_b, i_a)
        WRITE(16, '(F15.11)') q1_20(i_b, i_a)
        WRITE(114, '(F15.11)') h1_20(i_b, i_a)
        WRITE(115, '(F15.11)') gn1_20(i_b, i_a)                     
        WRITE(116, '(F15.11)') ct1_20(i_b, i_a)                     
    end do
    end do
    
    CLOSE(10)
    CLOSE(11)
    CLOSE(12)
    CLOSE(13)
    CLOSE(16)
    CLOSE(114)
    CLOSE(115)
    CLOSE(116)

    open (114, FILE='graphs/commit/period1/haut.txt',STATUS='UNKNOWN')
    open (115, FILE='graphs/commit/period1/gnaut.txt',STATUS='UNKNOWN')
    open (116, FILE='graphs/commit/period1/cttaut.txt',STATUS='UNKNOWN')
    
    do i_a = 1,na2 
        WRITE(114, '(F15.11)') haut1_20(i_a)
        WRITE(115, '(F15.11)') gnaut1_20(i_a)                     
        WRITE(116, '(F15.11)') ctaut1_20(i_a)                     
    end do
    
    CLOSE(114)
    CLOSE(115)
    CLOSE(116)

PRINT*, 'VALUES FOR COMMITMENT T=1 SOLVED WITH PERIOD-2 CUTS'

end subroutine val_COMMIT023
!**********************************************************************************
subroutine val_COMMIT022

USE GLOBAL
!USE OMP_LIB

IMPLICIT NONE

INTEGER :: d

REAL(8)::a_initial,b_initial
INTEGER::i_default_global
REAL(8)::a_valor,b_valor,b_valor2,v_valor,v_valor2,q_fun,b_next,b0_next
REAL(8)::b_next1(na2) 

REAL(8)::q,interpolate,b_next0
REAL(8)::FDATA(nb),DLEFT,DRIGHT,BREAK_GRID(nb),CSCOEF(4,nb),DCSVAL

INTEGER i_b, i_a, i_g,i, j, ILEFT, IRIGHT,NINTV,DNORIN
REAL(8)::x0,x1(na),b_fun,bmin_feas,new_value
INTEGER::index_bfeas(na),index_opt,ii_b,iii_b,nbaux,load_results
INTEGER::i_gmax,i_gmax0(1),i_a_index,i_b_index
CHARACTER(len=29)::filename
REAL(8),DIMENSION(nb,na)::q_matrix_nodef_2,v_matrix_2,q_scheme_2,v0_matrix_2,vaut_matrix_2,probdef_matrix_2
REAL(8),DIMENSION(nb2,na2)::v_matrix_20,q_scheme_20,q_matrix_nodef_20
INTEGER,DIMENSION(nb,na)::index_2
REAL(8)::objective_fun,objective_fun_2,interpolate3,b00,probdef_fun
EXTERNAL::objective_fun,objective_fun_2,probdef_fun

load_results = 1
PRINT*,'PERIOD 1: UPLOADING VALUE FCNS AND PRICES FOR PERIOD 2'

!upload variables for period 2
IF (load_results.EQ.1) THEN
    open (10, FILE='graphs/commit/period2/v_matrix.txt',STATUS='OLD')
    !open (15, FILE='graphs/commit/period2/probdef_matrix.txt', STATUS='OLD')
    open (16, FILE='graphs/commit/period2/q_paid_matrix.txt', STATUS='OLD')
      
    do i_b = 1,nb
    do i_a = 1,na                
        READ(10, '(F15.8, X, F15.8, X, F15.8)') v_matrix_2(i_b, i_a), v0_matrix_2(i_b, i_a), vaut_matrix_2(i_b, i_a)
        !READ(15, '(F15.11)') probdef_matrix_2(i_b, i_a)
        READ(16, '(F15.11, X, F15.11)') q_matrix_2(i_b, i_a), q_matrix_nodef_2(i_b, i_a)
    end do
    end do
    CLOSE(10)
    !CLOSE(15)
    CLOSE(16)
ENDIF
       
!to compute value fcns and prices correctly, set
v0_matrix(:,:)  =v0_matrix_2(:,:)
vaut_matrix(:,:)=vaut_matrix_2(:,:)
q_matrix_nodef(:,:) = q_matrix_nodef_2(:,:)
        
!to compute spline for feasible debt gridpoints
do i_a = 1,na
    new_value=-1d5
    i_b=0
    index_opt = 0
        
    do while (new_value <= -1d3)
        index_opt = i_b
        IF (i_b.EQ.nb) EXIT
        i_b=i_b+1
        new_value = v0_matrix(i_b,i_a)
    end do
    index_bfeas4(i_a)=index_opt     
ENDDO

index_bfeas2=MAXVAL(index_bfeas4)

do i_a = 1,na         
    IF (index_bfeas4(i_a)+1.LE.nb-1) THEN      !need at least to points for splines
        ii_b=index_bfeas4(i_a)+1
        nbaux =nb-ii_b+1
        FDATA(1:nbaux) = v0_matrix(ii_b:nb,i_a)
	    ILEFT  = 0
	    IRIGHT = 0
	    DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(ii_b+1) - bgrid(ii_b))
	    DRIGHT = (FDATA(nbaux)-FDATA(nbaux-1)) / (bgrid(nbaux) - bgrid(nbaux-1))
	    call  DCSDEC (nbaux, bgrid(ii_b:nb), FDATA(1:nbaux), ILEFT, DLEFT, IRIGHT, DRIGHT, break_matrix(i_a, ii_b:nb), coeff_matrix(i_a, :,ii_b:nb))
    ENDIF         
ENDDO
         
do i_a = 1,na
    FDATA = q_matrix_nodef(:,i_a)
	ILEFT  = 0
	IRIGHT = 0
    DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(2) - bgrid(1))
	DRIGHT = (FDATA(nb)-FDATA(nb-1)) / (bgrid(nb) - bgrid(nb-1))
	call  DCSDEC (nb, bgrid, FDATA, ILEFT, DLEFT, IRIGHT, DRIGHT, BREAK_GRID, CSCOEF)
    break_matrix_q(i_a, :) = BREAK_GRID
    coeff_matrix_q(i_a, :,:) = CSCOEF
end do

IF (psi1.LE.500d0) THEN
    !! Compute the value of repayment
    LOOP_A1: do i_a = 1,na2
        a_initial   = agrid2(i_a)
        b_initial   = bgrid2(1)
                  
        i_default_global=2 !Country defaults
        b_valor = 0d0

        !choose gN and transfers optimally in autarky
        v_valor= -objective_fun(b_valor,a_initial,b_initial,i_default_global)
        haut1_20(i_a)  = x_global
        gnaut1_20(i_a) = gn_global
        ctaut1_20(i_a) = ctrad_global
        b_next1(i_a)   = b_valor
        vaut1_20(i_a)  = v_valor
    end do LOOP_A1
ELSE
    haut1_20  = 0d0
    gnaut1_20 = 0d0
    ctaut1_20 = 0d0
    b_next1    = 0d0
    vaut1_20  = -1d10
ENDIF

i_a_index =index_aglobal
i_b_index = index_global 

PRINT*,'i_b_index',i_b_index

    !!!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(a_initial,b_initial,i_default_global,v_valor,b_next0)
    !!!$OMP DO COLLAPSE(2)
    LOOP_A2: do i_a = i_a_index,i_a_index 
    LOOP_B:  do i_b = i_b_index,i_b_index 
                  a_initial = agrid2(i_a)
                  b_initial = bgrid2(i_b)

                  i_default_global=1  !Country does not default

                  !slash spending and set transfers to equilibrium values in first period for first-period cuts
                  constant_gn = (1d0-sizecut)*LINEAR_INT(bgrid,gn_base(11,:),b_initial) 
                  transf_global = LINEAR_INT(bgrid,transf_base(11,:),b_initial) 
                  call optimize_gN2(b_next0, v_valor,a_initial,b_initial,i_default_global)
                  
                vcon1_20(i_b, i_a) = v_valor             
                x0 = x_global

                gn1_20(i_b, i_a) = gn_global   
                h1_20(i_b, i_a)  = x0                     
                ct1_20(i_b, i_a) = ctrad_global   
                b1_20(i_b, i_a) = b_next0        
                
                if (vaut1_20(i_a) > vcon1_20(i_b, i_a)) THEN !IT IS OPTIMAL TO DEFAULT
                    def1_20(i_b,i_a) = 0
                    v1_20(i_b, i_a) = vaut1_20(i_a)
                else
                    def1_20(i_b, i_a) = 1
                    v1_20(i_b, i_a) = vcon1_20(i_b, i_a)
                END if
                
                q1_20(i_b, i_a) = q_fun(b1_20(i_b, i_a),a_initial) 
    end do LOOP_B
    end do LOOP_A2
    !!!$OMP END DO
    !!!$OMP END PARALLEL 
        
    !save solution for first period
    open (10, FILE='graphs/commit/period1/v.txt',STATUS='UNKNOWN')
    open (11, FILE='graphs/commit/period1/default.txt',STATUS='UNKNOWN')
    open (12, FILE='graphs/commit/period1/bb_next.txt',STATUS='UNKNOWN')
    open (13, FILE='graphs/commit/period1/b_next.txt',STATUS='UNKNOWN')
    open (16, FILE='graphs/commit/period1/q_paid.txt', STATUS='UNKNOWN')
    open (114, FILE='graphs/commit/period1/x.txt',STATUS='UNKNOWN')
    open (115, FILE='graphs/commit/period1/gn.txt',STATUS='UNKNOWN')
    open (116, FILE='graphs/commit/period1/ctt.txt',STATUS='UNKNOWN')
    open (117, FILE='graphs/commit/period1/ign.txt',STATUS='UNKNOWN')
    
    do i_b = 1,nb2
    do i_a = 1,na2
        WRITE(10, '(F15.8, X, F15.8, X, F15.8)') MAX(-99d0,v1_20(i_b, i_a)), &
        MAX(-99d0,vcon1_20(i_b, i_a)), MAX(-99d0,vaut1_20(i_a))
        WRITE(11, '(I7)') def1_20(i_b, i_a)
        bp1_20(i_b,i_a)    = LOCATE2(bgrid2,b1_20(i_b, i_a))
        
        WRITE(12, '(F15.11)') b1_20(i_b, i_a)
        WRITE(13, '(I7)') bp1_20(i_b, i_a)
        WRITE(16, '(F15.11)') q1_20(i_b, i_a)
        WRITE(114, '(F15.11)') h1_20(i_b, i_a)
        WRITE(115, '(F15.11)') gn1_20(i_b, i_a)                     
        WRITE(116, '(F15.11)') ct1_20(i_b, i_a)                     
    end do
    end do
    
    CLOSE(10)
    CLOSE(11)
    CLOSE(12)
    CLOSE(13)
    CLOSE(16)
    CLOSE(114)
    CLOSE(115)
    CLOSE(116)

    open (114, FILE='graphs/commit/period1/haut.txt',STATUS='UNKNOWN')
    open (115, FILE='graphs/commit/period1/gnaut.txt',STATUS='UNKNOWN')
    open (116, FILE='graphs/commit/period1/cttaut.txt',STATUS='UNKNOWN')
    
    do i_a = 1,na2 
        WRITE(114, '(F15.11)') haut1_20(i_a)
        WRITE(115, '(F15.11)') gnaut1_20(i_a)                     
        WRITE(116, '(F15.11)') ctaut1_20(i_a)                     
    end do
    
    CLOSE(114)
    CLOSE(115)
    CLOSE(116)

PRINT*, 'VALUES FOR COMMITMENT T=0 SOLVED FOR CURRENT CUTS'

end subroutine val_COMMIT022
!*****************************************************************************
SUBROUTINE INITIAL_COMMIT20

USE GLOBAL

IMPLICIT NONE

REAL(8)::xval(6),fvalor,func
EXTERNAL::func

INTEGER::i,j,i_a,i_b,i0(1)
REAL(8)::FDATA(nb),DLEFT,DRIGHT,BREAK_GRID(nb),CSCOEF(4,nb),DCSVAL,epsilongrid(ncdf), dev_q_nodef,new_value
INTEGER::index_bfeas(na),index_opt,ii_b,nbaux
REAL(8)::q_matrix_nodef2(nb, na)
INTEGER indice_b, ILEFT, IRIGHT,NINTV,DNORIN,new_index_b,index_b !i_def_opt

!upload baseline in period 1-2
open (10, FILE='graphs/commit_base/period1/gvcon.txt')
open (11, FILE='graphs/commit_base/period1/gdef1.txt')
open (12, FILE='graphs/commit_base/period1/gbnext.txt')
open (13, FILE='graphs/commit_base/period1/gibp.txt')
open (16, FILE='graphs/commit_base/period1/gq_nodef.txt')
open (114, FILE='graphs/commit_base/period1/gh.txt')
open (115, FILE='graphs/commit_base/period1/ggn.txt')
open (116, FILE='graphs/commit_base/period1/gct.txt')
    
!make sure that indice order is correct
do i_a = 1,na2
do i_b = 1,nb2
    READ(10, *) vcon1(i_a,i_b)
    READ(11, '(I7)') def1(i_a,i_b)
    READ(12, *) bbp1(i_a,i_b)
    READ(13, '(I7)') bp1(i_a,i_b)
    READ(16, *) q1(i_a,i_b)                
    READ(114, *) h1(i_a,i_b)
    READ(115, *) gn1(i_a,i_b) 
    READ(116, *) ct1(i_a,i_b)                     
end do
end do
  
CLOSE(10)
CLOSE(11)
CLOSE(12)
CLOSE(13)
CLOSE(16)
CLOSE(114)
CLOSE(115)
CLOSE(116)

open (113, FILE='graphs/commit_base/period1/gvaut.txt')
open (114, FILE='graphs/commit_base/period1/ghaut.txt')
open (115, FILE='graphs/commit_base/period1/ggnaut.txt')
open (116, FILE='graphs/commit_base/period1/gctaut.txt')
    
do i_a = 1,na2 
    READ(113, *) vaut1(i_a)
    READ(114, *) haut1(i_a)
    READ(115, *) gnaut1(i_a)                     
    READ(116, *) ctaut1(i_a)                     
end do
    
CLOSE(113)
CLOSE(114)
CLOSE(115)
CLOSE(116)



PRINT*,'POLICIES FOR BASELINE UPLOADED'
!!***************************************************************
do i=1,na2
   agrid2(i) = amin + (amax - amin) * REAL((i-1d0) / (na2 - 1d0))
end do

CALL BUILD_TRANSITION(na2,agrid2,rhoa,sigmaa,mua,prob)
!***************************************************************** 
!upload solution in period 1 for spending cuts 
open (10, FILE='graphs/commit/period1/v.txt')
open (11, FILE='graphs/commit/period1/default.txt')
open (12, FILE='graphs/commit/period1/bb_next.txt')
open (13, FILE='graphs/commit/period1/b_next.txt')
open (16, FILE='graphs/commit/period1/q_paid.txt')
open (114, FILE='graphs/commit/period1/x.txt')
open (115, FILE='graphs/commit/period1/gn.txt')
open (116, FILE='graphs/commit/period1/ctt.txt')
open (117, FILE='graphs/commit/period1/ign.txt')
    
    do i_b = 1,nb2
    do i_a = 1,na2
        READ(10, '(F15.8, X, F15.8, X, F15.8)') v1_20(i_b, i_a),vcon1_20(i_b, i_a),vaut1_20(i_a)
        READ(11, '(I7)') def1_20(i_b, i_a)
        READ(12, '(F15.11)') b1_20(i_b, i_a)
        READ(13, '(I7)') bp1_20(i_b, i_a)
        READ(16, '(F15.11)') q1_20(i_b, i_a)   
        READ(114, '(F15.11)') h1_20(i_b, i_a)
        READ(115, '(F15.11)') gn1_20(i_b, i_a) 
        READ(116, '(F15.11)') ct1_20(i_b, i_a)                     
    end do
    end do

CLOSE(10)
CLOSE(11)
CLOSE(12)
CLOSE(13)
CLOSE(16)
CLOSE(114)
CLOSE(115)
CLOSE(116)

open (114, FILE='graphs/commit/period1/haut.txt')
open (115, FILE='graphs/commit/period1/gnaut.txt')
open (116, FILE='graphs/commit/period1/cttaut.txt')
    
do i_a = 1,na2 
    READ(114, '(F15.11)') haut1_20(i_a)
    READ(115, '(F15.11)') gnaut1_20(i_a)                     
    READ(116, '(F15.11)') ctaut1_20(i_a)                     
end do
    
CLOSE(114)
CLOSE(115)
CLOSE(116)

do i_b = 1,nb2
do i_a = 1,na2
    if (vaut1(i_a) > vcon1(i_a,i_b)) then
        def1(i_a,i_b) =0
    else
        def1(i_a,i_b) =1
    end if
enddo
enddo
             
open (10, FILE='graphs/commit_base/period1/v.txt')
open (16, FILE='graphs/commit/period2/q_paid_matrix.txt')
open (12, FILE='graphs/commit_base/period1/b_next.txt')
      do i_b = 1,nb
         do i_a = 1,na

             READ(10, '(F18.10, X, F18.10, X, F15.10)') v_matrix(i_b, i_a), &
                     v0_matrix(i_b, i_a), vaut_matrix(i_b, i_a)
            READ(16, '(F15.11, X, F15.11)') q_matrix_2(i_b, i_a), q_matrix_nodef(i_b, i_a)
            READ(12, '(F15.11)') b_next_matrix(i_b, i_a)
             if (vaut_matrix(i_b, i_a) > v0_matrix(i_b, i_a)) then
                 default_decision(i_b, i_a) =2
             else
                 default_decision(i_b, i_a) =1
             end if
         end do
      end do
           
CLOSE(12)
CLOSE(16)
CLOSE(17)

!to compute spline for feasible debt gridpoints
do i_a = 1,na
    new_value=-1d5
    i_b=0
    index_opt = 0
        
    do while (new_value <= -1d3)
        index_opt = i_b
        IF (i_b.EQ.nb) EXIT
        i_b=i_b+1
        new_value = v0_matrix(i_b,i_a)
    end do
    index_bfeas4(i_a)=index_opt     
ENDDO

index_bfeas2=MAXVAL(index_bfeas4)

do i_a = 1,na         
    IF (index_bfeas4(i_a)+1.LE.nb-1) THEN      !need at least to points for splines
        ii_b=index_bfeas4(i_a)+1
        nbaux =nb-ii_b+1
        FDATA(1:nbaux) = v0_matrix(ii_b:nb,i_a)
	    ILEFT  = 0
	    IRIGHT = 0
	    DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(ii_b+1) - bgrid(ii_b))
	    DRIGHT = (FDATA(nbaux)-FDATA(nbaux-1)) / (bgrid(nbaux) - bgrid(nbaux-1))
	    call  DCSDEC (nbaux, bgrid(ii_b:nb), FDATA(1:nbaux), ILEFT, DLEFT, IRIGHT, DRIGHT, break_matrix(i_a, ii_b:nb), coeff_matrix(i_a, :,ii_b:nb))
    ENDIF         
ENDDO
         
do i_a = 1,na
    FDATA = q_matrix_nodef(:,i_a)
	ILEFT  = 0
	IRIGHT = 0
    DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(2) - bgrid(1))
	DRIGHT = (FDATA(nb)-FDATA(nb-1)) / (bgrid(nb) - bgrid(nb-1))
	call  DCSDEC (nb, bgrid, FDATA, ILEFT, DLEFT, IRIGHT, DRIGHT, BREAK_GRID, CSCOEF)
    break_matrix_q(i_a, :) = BREAK_GRID
    coeff_matrix_q(i_a, :,:) = CSCOEF
end do

PRINT*,'POLICIES FOR SPENDING CUTS UPLOADED'

END SUBROUTINE INITIAL_COMMIT20
!!!*****************************************************************************